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A  method  is  described  for  assessing  the  performance  of  acoustic  arrays 
used  to  determine  source  bearings.  The  method  involves  calculation  of  the 
Cramer-Rao  lower  bound  (CRLB),  which  characterizes  the  best  performance 
obtainable  for  a  given  array  configuration  and  set  of  operating  conditions. 
The  CRLB  calculations  are  used  to  show  that  the  performance  of  the  arrays 
depends  on  the  sensor  configuration,  the  acoustic  frequency,  the  distance 
from  the  source,  the  background  noise,  and  atmospheric  turbulence.  Near 
to  the  source,  and  at  low  frequencies,  it  is  the  background  noise  that  limits 
array  performance.  Turbulence  becomes  the  limiting  factor  as  the  distance 
and  frequency  are  increased.  When  performance  is  limited  by  turbulence, 
the  calculations  are  found  to  be  very  sensitive  to  the  particular  turbulence 
model  used.  The  von  Karman  type  of  model  appears  to  provide  the  most 
reasonable  CRLB  calculations.  The  turbulent  degradation  is  found  to  result 
primarily  from  small-scale  fluctuations  in  the  wind  velocity.  Performance 
predictions  for  ARL’s  Remote  Netted  Acoustic  Detection  System  (RNADS) 
are  presented;  it  is  found  that  RNADS  should  provide  angle-of-arrival  (AO A) 
accuracy  of  several  degrees  or  better  for  signal-to-noise  ratios  above  10  dB 
during  most  atmospheric  turbulence  conditions. 
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1.  Introduction 


Acoustic  sensor  arrays  can  be  used  to  detect,  locate,  and  classify  targets 
in  atmospheric  battlespaces.  Their  main  advantages,  in  comparison  to  radar 
and  optical  systems,  are  low  cost,  functionality  in  non-line-of-sight  situations, 
and  insusceptibility  to  simple  countermeasures  (Srour  and  Robertson,  1995). 
These  advantages  help  explain  why  acoustical  systems  have  achieved  a  surge 
in  popularity  in  recent  years;  some  of  these  new  acoustical  systems  are  listed 
in  figure  1.  The  deployment  of  an  acoustical  array  for  target  tracking  is 
illustrated  in  figure  2.  The  operational  principle  is  to  use  the  phase  differences 
between  the  signals  recorded  at  the  individual  microphones  to  determine  the 
orientation  of  the  wavefronts,  which  corresponds  to  the  bearing  of  the  target. 

Acoustical  systems  are  not  without  shortcomings.  Most  importantly,  the 
performance  of  acoustical  arrays  depends  strongly  on  meteorological  condi¬ 
tions.  It  is  normally  enhanced  by  still,  nighttime  conditions.  On  the  other 
hand,  the  atmospheric  turbulence  characteristic  of  windy  conditions  or  a 
sunny  day  can  cause  significant  degradation.  This  is  because  array  beam- 
forming  depends  on  good  mutual  coherence  of  the  signals  received  by  the 
individual  sensors;  turbulence  causes  a  loss  of  such  coherence. 

Despite  the  recent  popularity  of  acoustical  arrays  and  the  importance 
of  atmospheric  turbulence  in  determining  the  performance  of  these  systems, 
so  far  little  effort  has  been  devoted  to  quantifying  the  turbulence  effects. 
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Use  acoustic  arrays  to  detect,  classify, 
estimate  bearing,  and  cue  optics. 


Use  multiple,  networked  arrays  to 
track  targets  by  bearing  estimation 
and  triangulation. 


Figure  1:  Partial  listing  of  current  generation  of  acoustical  systems  used  to 
track  sources  on  battlefield. 


Figure  2:  Deployment  of  acoustical  array  for  tracking  target.  Phase  infor¬ 
mation  at  microphones  is  used  to  determine  orientation  of  incoming  sound 
wavefronts.  Turbulence  causes  random  fluctuations  in  wavefiront  orientation. 
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In  order  to  place  the  previous  research  in  perspective,  it  is  helpful  to  break 
down  the  general  topic  of  turbulence  effects  on  acoustic  sensor  arrays  into 
two  subtopics:  first,  the  effect  of  turbulence  on  reducing  the  mutual  coher- 
ence  between  a  single  pair  of  microphones,  as  a  function  of  their  separation; 
and  second,  the  statistical  analysis  of  the  performance  of  a  full  array,  given 
knowledge  of  the  coherence. 

Some  of  the  earliest  work  on  coherence  is  summarized  by  Tatarskii  (1961, 
1971).  Tatarskii’s  theoretical  treatment  has  served  as  the  basis  for  most 
subsequent  work  throughout  the  aooustical  spectrum.  The  main  difficulty 
with  Tatarskii’s  theory  is  its  inapplicability  to  large-scale  turbulence  (i.e., 
eddy  sizes  larger  than  several  meters),  which  has  a  complicated  anisotropic 
and  inhomogeneous  structure.  Acoustical  systems  that  operate  in  frequency 
ranges  of  several  hundred  hertz  and  lower  are  strongly  affected  by  large-scale 
turbulence. 

Tatarskii  also  described  some  experimental  evidence  for  his  treatment  of 
coherence,  although  the  studies  he  considered  dealt  mostly  with  ultrasonic 
frequencies,  which  are  less  sensitive  to  large-scale  turbulence  structure.  Ex¬ 
perimental  studies  on  the  coherence  of  audible-range  sound  waves  have  been 
performed  during  recent  decades  by  Daigle  et  al  (1983)  and  Havelock  et 
al  (1995). 

The  second  topic  mentioned  above,  the  statistical  analysis  of  array  per¬ 
formance,  has  apparently  not  been  addressed  previously  for  a  turbulent  at¬ 
mosphere.  Song  and  Ritcey  (1996),  however,  do  consider  the  performance  of 
acoustic  arrays  in  a  randomly  fluctuating  ocean.  Their  analysis  applies  quite 
well  to  the  atmospheric  case,  so  long  as  the  statistical  model  used  for  the 
fluctuations  is  appropriate  to  the  atmosphere.  The  analysis  in  this  report  is 
derived  in  part  from  Song  and  Ritcey’s. 

Section  2  provides  a  statistical  formulation  of  array  performance.  A  ma¬ 
trix  representation  for  the  sensor  cross  correlations  is  employed  to  cast  the 
equations  in  a  particularly  simple  format.  Modeling  of  turbulence  effects 
on  the  signal  coherence  is  discussed  in  section  3.  Section  4  presents  ex¬ 
ample  assessments  of  array  performance  for  realistic  meteorological  condi¬ 
tions.  Among  the  array  configurations  considered  is  the  one  used  by  the 
Army  Research  Laboratory  (ARL)  Remote  Netted  Acoustic  Detection  Sys¬ 
tem  (RNADS)  (Srour  and  Robertson,  1995).  RNADS  is  fairly  representative 
of  the  acoustical  systems  listed  in  figure  1. 

It  will  become  clear  from  many  of  the  example  calculations  in  this  report 
that  accurate  assessments  of  acoustic  array  performance  require  careful  at¬ 
tention  to  modeling  of  atmospheric  turbulence.  Therefore,  I  include  some  of 
the  relevant  turbulence  modeling  issues  and  mathematics  in  an  appendix. 


2. 


General  Matrix  Formulation 


2.1  Angle  of  Arrival  versus  Bearing 

The  analysis  described  in  this  section  is  for  two-dimensional  sensor  arrays. 
Three-dimensional  arrays  can,  of  course,  also  be  built  and  analyzed;  we  gain 
little  physical  insight,  however,  by  introducing  this  additional  complication. 

A  two-dimensional  array  can  be  used  to  determine  the  bearing  of  a  source 
only  in  the  plane  of  the  array.  For  example,  if  the  array  is  horizontal,  it  can  be 
used  to  determine  only  the  azimuthal  bearing  of  the  source.  A  vertical  array 
would  determine  the  elevation  angle.  Most  of  the  analysis  in  this  section 
applies  equally  well  to  horizontal  and  vertical  arrays.  However,  there  are 
practical  reasons  why  most  acoustical  systems  designed  for  the  battlefield 
are  horizontal  arrays,  rather  than  vertical  ones.  The  most  important  of  these 
is  refraction,  as  illustrated  in  figure  3.  Refraction  by  atmospheric  wind  and 
temperature  gradients  causes  the  orientation  of  the  wavefronts  to  gradually 
change.  Eventually,  there  may  be  little  resemblance  between  the  angle  of 
arrival  of  the  wavefronts  and  the  actual  bearing  of  the  source. 

When  propagation  is  upwind  (or  the  temperature  increases  with  height), 
sound  is  refracted  downward.  As  a  result,  there  can  be  a  ‘‘multipath  effect”: 
i.e.,  multiple  wavefronts  and  angles  of  arrival.  When  propagation  is  down¬ 
wind  (or  the  temperature  decreases  with  height),  sound  is  refracted  upward. 
Most  of  the  sound  energy  reaching  the  receivers  is  normally  scattered  there 
by  turbulence.  As  a  result,  the  angle  of  arrival  has  more  to  do  with  the 
location  of  the  dominant  scattering  volume  than  the  actual  bearing  of  the 
source. 

Fortunately,  horizontal  refraction  is  minimal  in  comparison  to  vertical 
refraction.  Hence  we  need  not  worry  about  this  complication  if  we  are  in¬ 
terested  only  in  horizontal  arrays.  It  is  worth  pointing  out,  however,  that  if 


wind 


Figure  3:  Effect  of  refraction  on  angle  of  arrival  of  acoustic  wavefronts.  As 
result  of  refraction,  actual  bearing  of  source  can  differ  from  angle  of  arrival. 


3 


we  have  access  to  propagation  models  that  can  determine  vertical  refractive 
effects  sufficiently  well,  it  should  be  possible  to  build  vertical  arrays  that  ac¬ 
curately  track  the  elevation  of  targets.  In  the  remaining  part  of  this  section, 
I  circumvent  the  issue  of  refraction  by  referring  to  the  accuracy  of  angle- 
of-arrival  (AOA)  estimates  for  the  acoustic  wavefronts,  as  opposed  to  the 
actual  angle  of  bearing  (AOB)  of  the  source.  The  distinction  between  AOA 
and  AOB  is  of  little  practical  significance  for  horizontal  tracking  arrays,  but 
should  be  kept  in  mind  for  vertical  tracking  arrays. 


2.2  Sensor  Cross-Correlation  Analysis 

Consider  a  planar  array  of  N  sensors.  The  signal  at  each  of  the  sensors  is 
assumed  to  have  two  contributions:  one  that  has  propagated  through  the 
atmosphere  from  a  source  and  arrived  at  the  array  with  wavefront  normal 
angle  and  the  second  consisting  of  random  noise.  These  two  contributions 
are  indicated  by  t)  and  n(f),  respectively.  The  signals  are  represented  in 
complex  notation,  as  is  common  in  acoustics  and  wave  propagation  generally. 
(See,  for  example,  Kinsler  et  al  (1982)  or  sect.  6.9  in  Burdic  (1984).  The 
complex  part  of  the  signal  may  be  regarded  as  the  Hilbert  transform  of 
the  real  signal.)  Boldface  is  used  here  to  indicate  matrices;  s(^,t)  and  n(^) 
are  column  vectors  having  N  elements,  each  element  corresponding  to  the 
signal  received  by  one  of  the  sensors.  The  source  contributions  vary  in  time 
owing  to  the  effect  of  random  turbulent  fluctuations  on  the  propagation.* 
The  noise  contributions  may  be  from  the  wind  or  other  acoustic  sources  on 
the  battlefield. 

The  total  received  signal  is 

p{ip,t)  =  s{'ip,t)  +  n(t).  (2.1) 

The  cross  correlations  between  the  sensor  signals  can  be  written  compactly 
in  matrix  form  as 


Rpp  {rp)  =  (p(^,  t)p{ip,  t)) ,  (2.2) 

where  the  angle  brackets  indicate  the  ensemble  average,  and  the  tilde  is  the 
conjugate  transpose.  Assuming  that  the  source  signal  and  noise  are  uncor¬ 
related, 

Rpp(^)=Rs,(V')  +  Rnn(V’)-  (2-3) 

Assuming  furthermore  that  the  noise  at  the  sensors  is  mutually  uncorrelated 
and  equal  in  variance,  we  can  write 

R^i^P)  =  Iiss{^P)  +  (Tll,  (2.4) 

*Of  course,  the  source  signal  may  vary  in  time  for  other  reasons.  For  example,  its  bearing 
and  distance  may  change  as  a  result  of  source  motion,  or  its  spectral  content  may  change  as 
a  result  of  mechanical  operations  (e.g.,  a  tank  accelerating).  Although  such  complications 
can  be  significant,  I  do  not  consider  them  in  this  report  in  order  to  preserve  the  focus  on 
turbulence  effects. 


where  is  the  noise  variance.  It  is  convenient  to  normalize  the  signals  and 
noise  so  that  s{tp,t)  has  unit  variance.  Then  cr^  becomes  a  noise-to-signal 
variance  ratio.  The  signal-to-noise  ratio  (SNR)  in  decibels  is  defined  as 

SNR  = -10  log  (2.5) 


2.3  Wavefront  Phase  Delays 

Realistic  modeling  of  {xp)  is  rather  diflficult.  While  this  correlation  matrix 
depends  most  importantly  on  the  relative  sensor  positions,  it  also  depends 
on  several  other  factors  that  affect  the  propagation,  such  as  turbulence,  re¬ 
fraction,  and  absorption  of  sound  by  air.  I  consider  the  turbulence  effects  in 
section  3;  for  now,  I  formulate  Rjs  {xp)  assuming  perfect  plane-wave  propa¬ 
gation.  This  sort  of  formulation  is  valid  if  the  dimensions  of  the  array  are 
small  compared  to  its  distance  from  the  source. 

Assuming  perfect  plane-wave  propagation,  the  elements  of  R^s  (V’)  de¬ 
pend  only  on  the  relative  phase  delays  between  the  incidence  of  the  acoustic 
wavefronts  at  the  receivers.  The  situation  is  shown,  for  two  sensors,  in  fig¬ 
ure  4. 

In  the  figure,  one  sensor  is  at  the  origin,  while  the  second  is  at  (x,y). 
The  AOB  of  the  second  sensor  is  a  =  arctan  (y/x).  Suppose  there  are  plane 
waves  incident  on  the  sensors,  having  an  orientation  angle  (3.  Note  that  ^  = 
xp—ir/2,xp  being  the  angle  of  the  wavefront  normal.  The  distance  between  the 
receivers,  projected  onto  the  wavefront  normal,  is  dx  =  dsin  {0  —  a),  where 
=  x'^  +  y'^.  Therefore,  the  phase  delay  between  wavefront  incidence  at  the 
two  receivers  is  kdx  =  kd sin  {fS  —  a)  =  —kd cos  {xp  —  a),  where  A:  =  27r/A  is 
the  wavenumber.  The  signal  at  sensor  two  is  equal  to  the  signal  at  sensor 
one  times  the  phase  factor  exp  (ikdx). 

The  same  geometrical  arguments  hold  when  multiple  receivers  are  con¬ 
sidered.  The  phase  delay  between  sensors  m  and  n  is,  therefore, 

Smn  {fp)  =  exp  [-ikdmn  COS  {xp  -  Omn)]  ,  (2-6) 

where 

dmn  ~  (®”  Xrn)  +  (j/n  ”  J/m)  i 


Figure  4:  Planar  wavefronts  incident  on  pair  of  receivers. 
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and 


amn  =  arctan  [(j/„  -  ym)  /  {xn  -  Xm)]  ■ 

The  matrix  formed  from  the  phase  delays,  S(t/>),  is  generally  referred 
to  as  the  array  steering  matrix.  In  the  nonturbulent  case,  (^)  =  S  (^). 
I  show  in  section  3,  however,  that  this  simple  relationship  no  longer  holds 
when  turbulence  is  considered. 


2.4  Assessment  of  Array  Performance 

We  can  quantify  the  performance  of  the  sensor  array  by  calculating  the  mean 
square  error  (mse)  where  ip  is  the  estimated  AOA.  By  way  of 

an  introduction,  I  calculate  the  mse  for  a  simple  example:  an  array  of  two 
sensors  separated  by  a  distance  d. 

Suppose  there  are  two  sensors  on  the  y-axis  at  y  =  —d/2  and  y  =  +d/2. 
If  we  take  the  phase  at  the  sensor  at  y  =  +d/2  as  the  reference  phase, 
equation  (2.1)  becomes  for  this  case 


Pi 

1 

ni 

1 

ni 

. 

ikd  cos{%p— TV  /  2) 

+ 

g— ifcdsiny 

n2 

When  there  is  no  noise,  we  can  determine  'ip  exactly  from  the  phase  difference 
A<p  between  pi  and  P2' 

A(p  =  Zpi  —  Zp2  =  ^  {P1P2)  —  kdsinip, 


and  hence 


tp  =  arcsin 


(2.7) 


(The  asterisk  denotes  the  complex  conjugate.)  When  noise  is  present,  it  is 
natural  to  estimate  ip  by  averaging  several  samples  of  A4>  =  Z  (P1P2)  •  It  can 
be  shown  that  the  real  and  imaginary  parts  of  pip2  are 

3?  [P1P2]  =  cos  a  +  n2x  +  nix  cos  a  +  niy  sin  a  +  nix‘n2x  +  niyU^y, 

(2.8) 


^  \piP^  =  —  sin  a  —  n^y  +  niy  cos  a  —  nix  sin  a  —  nixn2y  +  niyn2x, 

(2.9) 

where  a  =  —kd  sin  ip,  Ujx  —  3?  [n^],  and  Ujy  =  Q  [uj].  Note  that  the  variance 
of  Ujx  and  njy  must  individually  be  0-^/2,  in  order  for  the  variance  of  nj  to 
be  equal  to  ^ 

It  is  not  possible  to  derive  closed-form  results  for  ip  from  equations  (2.8) 
and  (2.9)  in  general.  However,  if  we  suppose  that  cr^,  kd,  and  ip  are  all  much 
less  than  one  (i.e.,  low  noise,  small  separation  compared  to  the  wavelength, 
and  near  broadside  incidence),  we  have 

\P1P2]  c:^l  +  n2x  +  nix, 
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^  [P1P2]  n2y  +  niy. 

Hence  the  phase  angle  A0  is  approximately 


A<^ 


arctan 


/—a  —  n2y  + 

\  1  +  n2x  +  '^ix  ) 


(-a  -  n2y  +  niy)  (1  -  n2x  -  nix) 
—a  —  n2y  +  niy 
kd'tp  —  n2y  +  niy. 


We  find  the  estimated  AOA  ip  by  averaging  several  (M)  samples  of  A(p  and 
dividing  by  kd: 


>1  ivi 


M 


«lv 


—  n: 


2y 


m=l 


kd 


where  the  superscripts  m  indicate  statistically  independent  samples  of  the 
noise.  The  mse  can  now  be  calculated; 


{kdMy 

1 

{kdM)' 

1 


M 


EK-"S)' 

.m=l 

r2M  {al/2) 


M{kd) 


2- 


(2.10) 


For  most  array  geometries  and  processing  methods,  it  is  unfortunately 
impossible  to  calculate  the  mse  in  closed  form.  Fortunately,  it  is  possible  to 
perform  calculations  on  a  computer  that  give  useful  information  on  the  mse. 
A  particularly  useful  characterization  of  the  mse  is  provided  by  the  Cramer- 
Rao  theorem.  According  to  this  theorem,  the  minimum  mse  is  the  inverse  of 
the  Fisher  information  matrix  J  {ip)  (Scharf,  1991): 


(^(ptp-  >  J  ^  (V’)  •  (2-11) 

The  right-hand  side  of  equation  (2.11)  is  called  the  Cramer-Rao  lower  bound 
(CRLB).  In  our  case,  since  we  have  only  one  parameter  to  be  determined. 
Ip,  the  information  matrix  is  actually  a  scalar.  (For  multiple  sources,  the 
information  matrix  would  have  multiple  elements.)  Song  and  Ritcey  (1996) 
show,  for  signals  having  a  joint-Gaussian  probability  distribution,  that 

JW  =  Mtr(B;,'^R^>^).  (2.12) 

In  the  above,  tr()  is  the  trace  of  the  matrix,  and  M  is  the  number  of  sta¬ 
tistically  independent  samples  used  to  estimate  the  correlation  matrix  of 
p{ip,t). 

It  must  be  kept  in  mind  that  the  Cramer-Rao  theorem  gives  the  minimum 
possible  error.  Since  noise  and  turbulence  effects  on  the  propagation  are  ran¬ 
dom,  the  actual  error  is  also  random.  The  error  furthermore  depends  on  what 
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method  is  used  to  estimate  the  AOA.  Among  the  more  common  methods 
are  sum-difference  beamforming,  the  Bartlett  beamformer,  and  maximum- 
likelihood  estimation.  The  essential  task  to  be  accomplished  is  this:  given  a 
number  of  independent  samples  of  p(t/),  t)  collected  at  a  discrete  set  of  times 
tm,  and  a  theoretical  model  for  (ip),  we  desire  an  estimate  for  ip.  Most 
methods  involve  estimating  the  correlation  matrix  of  p{ip,  t)  from  the  actual 
measurements.  For  example,  the  following  matrix  may  be  computed: 

j  M 

m=l 

Given  enough  independent  samples,  we  would  expect  Cpp  to  converge  to 
Rpp.  The  reason  the  number  of  samples  M  appears  in  the  calculation  of  the 
Fisher  information  matrix  is  that  the  fidelity  of  the  estimate  Cpp  improves 
with  more  data,  so  that  the  estimate  of  the  AOA  is  improved. 

For  the  acoustic  array  performance  evaluations  given  in  this  report,  I 
provide  only  calculations  of  the  CRLB.  The  calculated  errors  are  therefore 
intrinsically  optimistic  assessments  of  actual  array  performance.  However, 
it  is  worth  pointing  out  that  Song  and  Ritcey  (1996),  in  their  calculations 
of  turbulence  effects  on  AOA  estimates,  were  able  to  achieve  performance 
very  near  the  CRLB  using  the  maximum-likelihood  method.  Therefore,  the 
CRLB  calculations  are  useful. 

In  most  of  the  remainder  of  this  report,  when  I  speak  of  the  CRLB,  I 
actually  mean  its  square  root.  Since  the  square  root  has  linear  dimensions 
(in  degrees),  it  is  more  intuitive  to  deal  with  than  the  actual  CRLB  (in  units 
of  degrees  squared).  The  square  root  of  the  CRLB  can  be  thought  of  as  a 
lower  bound  on  the  standard  deviation  of  the  measurement  error. 


2.4.1  Example  1:  Two- Element  Array  To  illustrate  the  CRLB 
analysis,  I  return  to  the  simple  example  considered  earlier:  that  of  two  sensors 
separated  by  a  distance  d.  In  matrix  form,  the  product  of  the  wavenumber 
and  sensor  spacings  is 


[^dmn]  — 


0  kd 
kd  0 


The  relative  angles  between  the  sensors  are 


0  7r/2 

-tt/2  0 


Hence,  from  equations  (2.6)  and  (2.4), 


(^)  ~ 


e 


—ikd  cos(V>+7r/2) 


g-ifecicos('^&-Tr/2) 


l  +  a^ 
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which,  using  elementary  trigonometric  identities,  we  can  simplify  to 


(V")  — 


^ikd  sin  “tp 


ikd  simp 


In  order  to  calculate  the  CRLB,  we  still  need  the  inverse  of  Rpp  (^),  as  well 
as  its  derivative  with  respect  to  V*  (see  eq  (2.12)).  The  inverse  of  the  matrix 
is  simply 

1  +  al  -e-^kdsm^  - 
_^^ikd sin  Ip  J  ’ 


■p  — 1 
^pp 


(ip)  = 


(2  +  (^1) 


and  the  derivative  is 

dKpp  _ 
dtp 


0 

ikd  cos 0 


Multiplying  these  matrices,  we  have 
idKpp  _  ikd  cos  Ip 


r; 


■pp 


dxp  (tI  (2  +  <72) 


-1 

(1  +  al) 


{1  +  al) 

1 


and  therefore 

/  _^5Rppp_i5Rpp\  _  2  {kd)  cos^  xp 
dtp  ^  dtp  <72(2  +  <r2)  • 

Hence  the  CRLB  (eq  (2.11))  is 

j-l  (^)  =  ^n(2  +  <),_ 

^  ^  2M{kdfcos^tP 


(2.13) 


Clearly  the  two-element  array  performs  best  when  the  source  is  perpendic¬ 
ular  to  the  array  axis  (^  =  0).  In  fact,  when  the  source  is  in  a  line  with  the 
receivers  (•^  =  7r/2),  the  CRLB  becomes  singular.  Note  that  equation  (2.13) 
is  equivalent  to  equation  (2.10)  if  one  makes  the  approximations  <72  <  1 
and  ^  -C  1.  Also,  the  two-element  array  has  a  ±7r  ambiguity  in  the  AOA 
estimation  (e.g.,  it  has  no  intrinsic  way  to  distinguish  tp  =  0  from  tp  =  it). 
The  CRLB  calculation  actually  assumes,  implicitly,  that  this  information  is 
somehow  available. 

A  polar  plot  of  the  CRLB  is  presented  in  figure  5  (solid  line),  for  cr^  =  Q.l 
and  d  =  A/2  {kd  =  tt).  The  CRLB  appears  as  two  parallel  lines  when  plotted 
in  this  manner. 


2.4.2  Example  2:  Triangular  Array  A  more  complicated  case  is 
that  of  an  equilateral  triangular  array.  With  vertices  at  {0,d/y/Z),  {—d/2, 
-d/2y/Z),  and  (d/2, -d/2\/3),  the  sensor  spacing  matrix  is 


[fc<i7nn]  — 


0  kd  kd 
kd  0  kd 
kd  kd  0 


and  the  angular  orientation  matrix  is 


[<^mn]  — 


0  — 27r/3  — 7r/3 

7r/3  0  TT 

27r/3  TT  0 
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Hence, 


\  ^  g~2fcrfcos(t/?H-27r/3)  ^—ikdcos('4)-\-‘n/2>) 

gtfcdcos(t/^+27r/3)  j  ^  ^2  ^-ikdcos{ip) 

gt/:(icos(^+7r/3)  ^ikdcos{ip)  ^  ^2 

Analytical  calculation  of  the  CRLB  becomes  extremely  complicated  for  this 
case.  It  is  rather  simple,  though,  to  conmute  numerical  results  with  a  matrix 
computation  package  such  as  Matlab®.  Results  for  d  =  A/2  are  shown 
(dotted  line)  superimposed  on  the  polar  plot  for  the  two-element  line  array 
(fig-  5). 

The  CRLB  for  an  equilateral  triangle  is  a  circle,  meaning  that  there  are 
no  preferred  directions  for  detection.  This  is  somewhat  surprising:  since  the 
line  array  has  two  preferred  directions,  one  might  expect  the  triangular  array 
to  have  three  preferred  directions.  But  all  test  cases  that  I  have  computed 
for  regular  polygons  (such  as  a  hexagon  with  side  length  d  =  A/2,  results  for 
which  are  also  shown  in  fig.  5,  double-dash-dotted  line)  have  yielded  an  omni¬ 
directional  CRLB.  The  underlying  mathematical  reasons  for  this  are  unclear 
to  me.  The  CRLB  calculation  for  a  random  configuration  of  three  or  more 
sensors,  also  shown  in  figure  5  (dotted  line),  yields  an  ellipsoidal  response. 
(The  ellipsoid  depicted  is  a  typical  result;  other  random  configurations  yield 
ellipsoids  having  different  orientations  and  axes.) 

In  connection  with  figure  5,  it  should  be  stressed  that  the  calculations 
are  for  the  CRLB,  and  not  for  actual  phased-array  beam  patterns  or  for 
specific  signal-processing  algorithms.  As  discussed  in  section  2.2.4,  the  CRLB 
represents  the  theoretically  best  attainable  performance.  Only  with  good 
algorithms  can  one  actually  obtain  similar  results. 


(^)  ~ 


CRLB  transverse  to  array  (°) 


Figure  5:  Polar  plot  of  Cramer-Rao  lower  bounds  for  AOA  estimations. 
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3.  Incorporation  of  Turbulence  Effects 

Turbulence  causes  random  fluctuations,  in  both  phase  and  amplitude,  of 
the  received  signals.  These  random  fluctuations  would  not  interfere  with 
AOA  estimation  if  they  caused  equal  and  synchronized  fluctuations  in  phase 
and  amplitude  at  all  sensors.  The  problem  is  that  the  phase  and  amplitude 
fluctuations  at  the  individual  array  elements  are  unsynchronized  with  one 
another.  In  this  section,  I  explain  how  the  effect  of  this  imperfect  coherence 
on  the  sensor  array  can  be  modeled. 


3.1  Effective  Index-of-Refraction  Fluctuation 

Turbulent  fluctuations  in  air  temperature,  moisture,  and  wind  velocity  all 
affect  acoustic  propagation.  Both  the  temperature  and  moisture  fluctuations 
cause  the  acoustic  index  of  refraction  to  fluctuate.  Let  us  define  the  squared 
index  of  refraction  as 


e  =  —  1. 


(3.1) 


In  the  above,  n  =  Co/c  is  the  index  of  refraction,  c  being  the  actual  soimd 
speed  and  cq  its  ambient  value.  The  speed  of  sound,  in  terms  of  the  absolute 
temperature  T  and  specific  humidity  q,  is 

c  ~  20.02^r(l  +  0.62g). 

Normally,  the  fluctuations  in  the  sound  speed,  cf  =  c— cq,  are  small  compared 
to  the  ambient  value.  Therefore, 


n 


2 


^0 


4 


(co  +  dy  4  +  2^00^ 


The  squared  index  of  refraction  is  then 


e  ~ 


Co  ' 


(3.2) 


The  effect  of  velocity  fluctuations  is  similar  to  the  effect  of  sound  speed 
fluctuations.  If  we  make  a  parabolic  approximation  to  the  acoustic  wave 
equation,  it  can  be  shown  that,  in  the  presence  of  turbulent  velocity  fluctu¬ 
ations,  e  can  be  replaced  by  the  following  effective  value  (Ostashev,  1994): 

Ce// =  e-2tii/co,  (3.3) 


in  which  ui  is  the  fluctuation  of  the  velocity  component  parallel  to  the  di¬ 
rection  of  propagation. 

The  statistical  spatial  structure  of  the  effective  index  of  refraction,  and 
its  effects  on  acoustic  propagation,  is  considered  in  section  3.3. 
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3.2  Mutual  Coherence  Function 


The  main  effect  of  turbulence  on  tracking  arrays  is  to  reduce  the  coherence 
of  the  signal  received  between  the  individual  array  elements.  This  effect  is 
normally  quantified  by  the  transverse  mutual  coherence  function  (MCF).  If 
Pm{i)  and  Pn(0  are  the  signals  at  two  sensors  whose  separation  p  is  parallel  to 
the  wavefronts  (transverse  to  the  propagation  direction),  the  MCF  is  defined 
as 

p  ,  s  ^  (Pm(*)Pn(0) 

(Pm(t)pJk(i))' 

(It  is  assumed  that  the  average  mean  intensity  at  the  two  receivers  matches, 
i.e.,  (Pm(<)PmW)  =  {Pn{t)Pn{t))-) 

Assuming  the  turbulence  is  homogeneous  and  isotropic,  the  MCF  is  given 
by  (Tatarskii,  1971) 

F  (p)  =  exp  |_![^  (0)  -  hejf  (p)]| ,  (3.4) 

where  x  is  the  propagation  distance,  p^  =  {ym  —  2/n)^  +  (^m  —  Znf  ■,  and 
(p)  is  called  the  effective  two-dimensional  correlation  function  (TCF).*^ 
The  TCF  is  a  special  type  of  correlation  function  involving  the  effective  index 
of  refraction.  (See  the  definition  and  note  on  terminology  in  the  appendix.)  It 
has  two  terms,  one  of  which  results  from  the  actual  fluctuations  in  the  index 
of  refraction,  and  the  second  from  velocity  fluctuations  (Ostashev,  1994): 

beff  (p)  =  b{p)-\-  bn  (p)  •  (3.5) 


3.3  Turbulence  Models 


The  problem  of  modeling  the  turbulence  effect  on  the  sensors  has  now  been 
reduced  to  one  of  modeling  the  TCFs  of  the  turbulence  field,  b  (p)  and  6ii  (p). 
The  appendix  shows  how  equations  for  the  TCFs  can  be  developed  for  two 
particular  turbulence  models,  the  Gaussian  and  von  Karman.  Both  the  Gaus¬ 
sian  and  von  Karman  models  come  in  two  forms,  one  for  scalars  (e.g.,  the 


*  Instead  of  equation  (3.4),  the  formula  for  the  MCF  is  often  given  as  (Rytov,  Kravtzov 
and  Tatarskii,  1989) 

r(p)  =  exp J  [1  ~  Jo  («±/!>)]  ^  (0,/C2,K3) I , 

where  Jo  is  the  Bessel  function  of  the  first  kind,  and  /ci  =  «!  +  Using  the  identity 


Jo{z) 


1 

2^  Jo 


exp  {iz  cos  6)  d9^ 


and  the  definition  for  the  TCF  (see  appendix,  eq  (A-4)),  we  can  show  the  two  formulas  to 
be  equivalent. 

^ Since  the  array  configurations  in  this  text  2ire  all  horizontal,  we  could  set  =  2n 
and  p  —  ym  ”  Vn  at  this  point.  It  is  better  not  to  do  so,  however.  The  main  reason  is 
that  the  turbulence  modeling  (see  the  appendix)  must  be  done  in  three  dimensions  in 
order  to  be  valid;  furthermore,  the  results  in  this  section  are  not  simplified  when  written 
in  two-dimensional  form. 


actual  index  of  refraction),  and  one  for  vectors  (e.g.,  the  velocity  fluctua¬ 
tions).  Each  form  contains  two  parameters,  a  variance  and  a  length  scale. 

The  variance  parameter  for  the  scalar  forms,  in  either  the  Gaussian  or 
von  Karman  models,  is  seen  from  equation  (3.2)  to  be 

^2  =  /,2\  ^  (3.6) 

For  the  vector  case,  we  have  directly  from  equation  (3.1) 

=  iM.  (3.7) 

The  length  scale  parameters  used  in  the  Gaussian  and  von  Karman  mod¬ 
els  are  not  equivalent.  Therefore,  a  common  standard  for  length  scales  is 
needed  for  comparing  the  models.  The  standard  I  adopt  is  the  so-called 
integral  length  scale,  which,  as  discussed  in  the  appendix,  we  can  find  by 
integrating  the  standard  correlation  function  from  0  to  oo,  and  then  divid¬ 
ing  by  the  variance.  For  scalar  fields,  this  procedure  yields  a  single  integral 
length  scale  £.  Two  scales  can  be  defined,  however,  for  vector  fields:  one 
for  which  the  direction  of  integration  is  parallel  to  the  orientation  of  the 
velocity  components,  and  one  in  which  the  integration  is  perpendicular  to 
the  velocity  components.  I  call  these  two  cases  the  parallel  and  perpendicular 
integral  length  scales,  denoted  by  C\\  and  C^,  respectively.  In  a  homogeneous, 
isotropic  turbulent  field,  £||  =  2C±  (Batchelor,  1953). 

For  the  scalar  form  of  the  Gaussian  model,  the  equation  for  the  TCF, 
derived  in  the  appendix,  is 

with  the  length  scale  parameter  given  by 

L=^C.  (3.9) 

For  the  vector  form  of  the  Gaussian  model,  the  TCF  is  (see  appendix) 


with 


(3.11) 


The  main  advantage  of  the  Gaussian  model  is  its  analytical  convenience. 
Unfortunately,  it  is  rather  unrealistic,  because  most  of  its  energy  occurs  con¬ 
centrated  at  the  scale  L.  In  actuality,  turbulent  energy  spans  a  broad  range 
of  spatial  scales.  The  von  Karman  model  captures  this  characteristic  of  the 
turbulence  much  more  satisfactorily.  In  particular,  the  von  Karman  model 
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includes  a  realistic  inertial  subrange,  in  which  the  energy  decays  according 
to  Kolmogorov’s  —5/3  power  law  (1941).  For  the  scalar  form  of  the  von 
Karman  model,  the  equation  for  the  TCF,  derived  in  the  appendix,  is 


b{p)  = 


2y^aH 

3v^r  {u) 


(3.12) 


where  is  the  modified  Bessel  function  of  the  second  kind,  and  the  length 
scale  parameter  is 


e  = 


3r  ju) 

2v^r(i/  +  i/2)  * 


(3.13) 


The  result  for  vectors  is 


where 


6(p)  = 


21/6^2^ 

\/ir(i/) 


e  = 


r(i/) 


r>Cil  = 


2T{iy) 


V^r  (j/ +  1/2)  y/nr{iy  +  l/2) 


Cl- 


(3.14) 


(3.15) 


TCFs  for  the  two  turbulence  models  are  compared  in  figure  6.  The  re¬ 
sulting  MCFs  are  shown  in  figure  7.  One  notable  feature  of  the  curves  is 
that  the  Gaussian  and  von  Karman  models  are  fairly  similar.  The  Gaussian 
curves  for  the  MCF  decay  somewhat  more  quickly,  because  that  model  has 
less  small-scale  turbulent  kinetic  energy  than  the  von  Karman  model.  The 
fact  that  the  two  models  yield  similar  results  for  the  TCF  and  MCF,  despite 
the  more  realistic  description  of  the  inertial  subrange  turbulence  intrinsic 
to  the  von  Karman  model,  indicates  that  the  inertial  subrange  turbulence 
is  of  secondary  importance  in  this  particular  application.  Although  Gaus¬ 
sian  models  often  do  not  work  well  in  wave  propagation  applications,  in  this 
particular  application  their  usage  may  be  justifiable.  I  consider  this  matter 
again  in  section  4.1. 

Figure  8  shows  the  MCF  for  several  values  of  k'^xC:  1,  10,  and  100. 
This  plot  makes  it  clear  that  signal  coherence  is  reduced  dramatically  with 
increasing  range  and  decreasing  wavelength. 


3.4  Modification  of  Signal  Correlation  Matrix 

The  MCF  describes  the  coherence  between  two  sensors  when  the  incident 
wavefronts  are  parallel  to  the  sensor  axis.  Unfortunately,  there  is  no  similar 
result  available  for  arbitrary  orientation  between  the  wavefronts  and  sensors. 
I  assume  here  that  the  coherence  between  a  pair  of  sensors  in  such  cases  can 
be  estimated  by  multiplication  of  the  phase  delay,  which  would  occur  in 
the  absence  of  turbulence  by  the  MCF  for  parallel  incidence.  Hence  we  can 
include  the  turbulence  effect  by  setting 

Rss  W  =  S  (^)  ©  T, 
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Figure  6:  Comparison  of  transverse  correlation  functions  for  different  turbu¬ 
lence  models. 


Figure  7:  Comparison  of  mutual  coherence  functions  for  several  different 
turbulence  models.  Value  used  for  k^xC  was  10. 


Figure  8;  Mutual  coherence  function  for  several  values  of  k^xC.  Scalar  form 
of  Gaussian  spectrum  was  used. 
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where  S  is  the  steering  matrix  (eq  (2.6)), 

Iran  —  T  {dnui)  , 

and  the  symbol  ©  indicates  the  Hadamard  matrix  product  (simple  element- 
by-element  multiplication). 


4.  Example  Calculations  and  Discussion 

4.1  Long  Baseline — Sensor  Decorrelation  Trade-Off 

If  it  were  not  for  cost  considerations  and  turbulence  effects,  it  would  be  de¬ 
sirable  to  make  acoustic  tracking  arrays  as  large  as  possible.  This  is  evident, 
for  example,  from  the  equation  for  the  CRLB  of  a  two-element  array  (eq 
(2.13)).  The  CRLB  in  this  case  is  inversely  proportional  to  the  square  of  the 
sensor  spacing.  Increasing  the  array  dimensions  is  desirable  because  it  in¬ 
creases  the  phase  differences  between  the  sensors,  making  the  measurement 
less  sensitive  to  noise. 

On  the  other  hand,  large  arrays  are  undesirable  when  turbulence  is  con¬ 
sidered.  Generally,  one  wants  the  dimensions  of  the  array  small  compared  to 
the  size  of  the  most  energetic  turbulent  eddies,  so  that  there  is  good  mutual 
coherence  across  the  array. 

The  trade-off  between  array  baseline  and  turbulence  degradation  is  il¬ 
lustrated  by  the  CRLB  calculations  shown  in  figure  9.  The  calculations  are 
for  parallel  incidence  {'tp  =  0)  on  a  five-element  line  array,  at  a  normalized 
distance  kx  =  1000  from  the  source.  The  spacing  kd  (=  27rd/A)  and  strength 
of  the  turbulence  are  varied  in  the  figures.  The  normalized  integral  length 
scale  is  kC  —  IOtt,  and  the  noise-to-signal  variance  is  =  0.1  (SNR  =  10). 
A  scalar  von  Karman  turbulence  spectrum  was  used. 

For  weak  turbulence  {o^  ^  10“®),  array  performance  essentially  can  be 
improved  without  bound  by  an  increase  in  the  sensor  spacing.  Performance  is 
poor  for  small  sensor  spacings,  because  phase  difference  information  becomes 
lost  in  the  measurement  noise.  For  strong  turbulence  (cr^  ^  10“'^),  array  per¬ 
formance  is  degraded,  particularly  if  the  size  of  the  array  is  comparable  to 
the  integral  length  scale  of  the  turbulence.  When  this  occurs,  an  “error  reso¬ 
nance”  in  the  CRLB  is  evident.  If  array  dimensions  are  either  much  smaller 
or  much  larger  than  the  turbulence  length  scale,  however,  array  performance 
suffers  little. 

Figure  10  is  similar  to  9,  except  that  the  noise  variance  has  been  decreased 
to  cr^  =  0.01.  Because  there  is  less  noise  than  in  the  previous  example,  smaller 
arrays  provide  better  performance.  Although  the  error  resonance  for  poor 
array  performance  is  still  clearly  evident  on  the  figure,  now  some  turbulent 
degradation  of  array  performance  is  observed,  even  for  the  smaller  arrays. 

A  scalar  Gaussian  spectrum  for  =  0.01  is  shown  in  figure  11.  The 
error  resonance  in  strong  turbulence  is  sharper  in  the  Gaussian  spectrum, 
while  the  error  for  small  arrays  (2d/A  <  1)  is  somewhat  less.  This  difference 
results  from  a  deficiency  of  the  Gaussian  model,  discussed  earlier:  this  model 
does  not  realistically  capture  the  energy  in  the  small-scale  turbulence. 


4.2  Dependence  of  Array  Performance  on  Meteorolog¬ 
ical  Conditions 

In  order  to  estimate  lower  bounds  on  array  performance  for  actual  scenarios, 
we  need  values  of  the  turbulence  spectral  parameters  and  £.  Obtaining 
these  is  actually  quite  a  challenging  problem,  since  these  parameters  depend 
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Figure  9:  Effect  of  changing  sensor  spacing  on  AOA  estimation  errors:  mod¬ 
erate  noise,  von  Karman  turbulence  model  case.  Array  is  five-element  line 
array,  and  calculations  are  for  cr^  =  0.1,  kC  =  IOtt. 
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Figure  10:  Effect  of  changing  sensor  spacing  on  AOA  estimation  errors:  low- 
noise,  von  Karman  turbulence  model.  Same  as  figure  9,  except  that  noise 
variance  has  been  reduced  to  =  0.01. 
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Figure  11:  Effect  of  changing  sensor  spacing  on  AOA  estimation  errors:  low- 
noise,  Gaussian  turbulence  model.  Except  for  change  in  turbulence  model, 
this  is  same  as  figure  10. 


strongly  on  meteorological  and  topographical  conditions.  There  is  no  en¬ 
tirely  satisfactory  turbulence  model  available.  Wilson  and  Thomson  (1994) 
have  proposed  a  Gaussian  turbulence  model  for  acoustical  applications  that 
accounts  for  many  of  the  important  features  of  the  turbulence,  although, 
like  other  Gaussian  models,  it  is  unrealistic  for  the  small-scale  turbulence 
structure.  Since  no  comparable  von  Karman  model  is  available,  I  adopt  the 
Wilson  and  Thomson  model  for  the  remaining  calculations  in  this  report. 

Because  the  full  Wilson  and  Thomson  model  is  rather  complicated,  I 
provide  here  only  a  simplified  version.  The  model  has  three  independent 
physical  parameters: 

1.  The  surface  wind  stress  t.  The  wind  stress  increases  with  increasing 
wind  speed. 

2.  The  surface  temperature  flux  Q,  from  the  ground  to  the  air.  The  tem¬ 
perature  flux  is  large  when  the  ground  is  warmer  than  the  overlying 
air,  such  as  typically  occurs  when  the  sun  heats  the  ground. 

3.  The  boundary-layer  inversion  height  zi.  This  is  the  height  where  the 
transition  occurs  between  the  relatively  cool,  atmospheric  boundary 
layer  air,  and  warmer  air  aloft  (the  free  troposphere).  In  fair  weather 
conditions,  this  boundary  is  quite  frequently  marked  by  stratocumulus 
(pufiy,  flattened)  clouds  at  several  hundred  to  several  thousand  meters 
altitude. 

It  is  convenient  to  define  several  scaling  parameters  in  terms  of  the  phys¬ 
ical  ones,  T,  Q,  and  Zi.  The  friction  velocity  is  defined  as 

(4.1) 

where  p  is  the  air  density.  The  surface-layer  temperature  scale  is 

r,  =  -— .  (4.2) 

u* 

The  free-convection  velocity  scale  is 

w,  =  (Qgzi/T,)^^^  =  {-u,ngzi/Ts)^/^ ,  (4.3) 

where  g  is  gravitational  acceleration  and  Tg  the  surface  temperature. 

The  Wilson  and  Thomson  model  includes  four  additive  contributions: 
small-scale  temperature  fluctuations,  small-scale  velocity  fluctuations,  tem¬ 
perature/velocity  covariance,  and  large-scale  velocity  fluctuations.  The  co- 
variance  term  is  neglected  in  this  report.  The  small-scale  temperature  (scalar) 
variance  is  modeled  as 


The  small-scale  velocity  (vector)  variance  is 
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For  both  of  the  small-scale  contributions,  the  corresponding  length  scales 
simply  equal  the  height: 


L  =  z. 

The  large-scale  velocity  (vector)  variance  is 


with  corresponding  length  scale 


L  =  OAzi. 

Prom  equation  (3.8),  the  resulting  scalar  TCF  is 

while  from  equation  (3.10)  the  vector  TCF  is 

6„(p)  =  8.4z  (^)'«<p  (-g]  +0.14Z,  (^) 


exp  - 


O.Wzf 


(4.4) 


(4.5) 


The  effective  TCF  is  simply  the  sum  of  b{p)  and  bn{p)  (see  eq  (3.5)). 
It  is  plotted  in  figure  12  for  u*  =  0.5  m/s,  T*  =  —0.3  K,  z  =  2  m,  and 
Zi  =  1000  m.  These  parameters  are  representative  of  a  very  sunny,  some¬ 
what  windy  afternoon,  what  atmospheric  scientists  would  call  a  convective 
boundary  layer. 


Figure  12:  Two-dimensional  correlation  function  (TCF)  for  convective,  at¬ 
mospheric  boundary  layer  conditions  (u*  =  0.5  m/s,  T*  =  —0.3  K,  z  =  2  m, 
and  Zi  =  1000  m). 


The  figure  clearly  shows  that  the  large-scale  velocity  fiuctuations  (the 
second  term  in  eq  (4.5))  dominate  the  effective  TCF.  The  reason  the  TCF 
is  dominated  by  large-scale  structure  is  evident  from  the  TCF  definition  (eq 
(A-4)):  by  setting  ki  =  0  before  performing  the  two-dimensional  (2D)  inverse 
Fourier  transform,  we  filter  out  the  small-scale  structure. 

The  observation  that  large-scale  structure  dominates  the  TCF  can  lead 
to  the  incorrect  conclusion  that  large-scale  turbulence  is  most  important  for 
determining  the  MCF.  Calculations  of  the  MCF  (eq  (3.4))  actually  depend 
on  the  difference  beff  (0)  —  be/f  (p),  not  simply  6e//  (p))  ^  was  plotted  in 
figure  12.  This  difference  function  is  plotted  for  the  convective  boundary  layer 
case  in  figure  13.  Now  we  see  that,  in  actuality,  the  small-scale  wind  velocity 
fluctuations  are  most  significant  for  MCF  calculations  when  p  <  10  m,  the 
situation  of  primary  interest  for  acoustical  arrays.  The  large-scale  velocity 
fluctuations  become  most  important  only  when  p  >  100  m.  Temperature 
fiuctuations  are  unimportant  compared  to  velocity  at  all  separations. 

The  difference  bgjf  (0)— 6e//  (p)  is  plotted  for  two  other  cases  in  figures  14 
and  15.  Figure  14  is  for  =  0.5  m/s  and  T*  =  —0.01  K,  with  z  and  zi  the 
same  as  earlier.  These  parameters  are  representative  of  what  atmospheric  sci¬ 
entists  call  neutral  conditions,  such  as  normally  occur  on  an  overcast,  windy 
day.  Figure  15  is  for  u*  =  0.1  m/s  and  r»  =  0.01  K,  representative  of  calm 
conditions.  In  these  two  cases,  as  before,  small-scale  velocity  fluctuations  are 
dominant  when  p  <  10. 

Example  calculations  of  the  MCF  as  a  function  of  u*  and  T*  are  shown  in 
figure  16.  The  calculations  are  for  p  equal  to  one-half  wavelength,  /  =  200  Hz, 
z  =  2  m,  Zi  —  1000  m,  and  propagation  distance  x  =  500  m.  The  resulting 
CRLB  is  shown  in  figure  17  for  a  five-element  line  array,  with  spacing  A/2, 
at  parallel  incidence  (■)/)  =  0).  The  number  of  samples  was  chosen  as  M  =  5, 
and  the  noise  variance  was  low,  cr^  =  0.01  (SNR  =  20). 


Figure  13;  TCF  difference  function  b^ff  (0)-6e//  (p)  for  convective  boundary 
layer  conditions. 
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Figure  14:  TCF  difference  function  be/f  (0)  -  beff  {p)  for  nearly  neutral  con¬ 
ditions:  w*  =  0.5  m/s,  r*  =  -0.01  K,  2:  =  2  m,  and  Zi  =  1000  m. 
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Figure  15:  TCF  difference  function  beff  {0)  —  beff  (p)  for  calm  boundary  layer 
conditions:  tt*  =  0.1  m/s,  T*  =  0.01  K,  ^  —  2  m,  and  zt  =  1000  m. 


Figure  16:  MCF  calculations  as  a  function  of  turbulence  conditions.  Propa¬ 
gation  distance  500  m,  frequency  200  Hz,  and  spacing  one-half  wavelength. 


Figure  17:  CRLB  as  a  function  of  turbulence  conditions,  for  a  five-element 
line  array  at  parallel  incidence.  Propagation  distance  500  m,  frequency  200 
Hz,  and  sensor  spacing  one-half  wavelength. 


4.3  Calculations  for  Circular  RNADS  Array 

The  purpose  of  ARL’s  RNADS  is  to  acoustically  detect  ground  and  air  tar¬ 
gets  in  battlefield  environments  (Srour  and  Robertson,  1995).  The  normal 
configuration  for  the  acoustic  array,  shown  in  figure  18,  has  seven  sensors: 
six  spaced  equally  along  a  circle,  and  a  seventh  at  the  origin.  In  this  sec¬ 
tion,  I  examine  the  theoretically  best  obtainable  performance  (CRLB)  by 
the  RNADS  array  in  various  turbulent  environments.  CRLB  calculations 
are  provided  for  the  same  three  representative  cases  of  turbulence  parame¬ 
ters  discussed  in  section  4.2:  convective  («*  =  0.5  m/s,  T*  =  —0.3  K),  neutral 
(n,  =  0.5  m/s,  T*  =  —0.01  K),  and  calm  (tt*  =  0.1  m/s,  T.  =  —0.01  K). 
In  all  cases  1  set  z  =  2  m,  Zi  =  1000  m,  and  M  =  5.  Note,  as  discussed 
in  section  2,  that  regular  polygons  such  as  the  RNADS  array  have  a  CRLB 
that  is  independent  of  the  AOA  V’-  Therefore,  it  is  unnecessary  to  plot  the 
CRLB  as  a  function  of  ^  for  RNADS. 

Two  sets  of  calculations  for  convective  conditions  are  shown  in  figures  19 
and  20,  as  a  function  of  the  propagation  distance  and  the  acoustic  frequency. 
The  difference  between  the  figures  is  the  value  used  for  the  noise  variance: 

=  0.01  (SNR  =  20  dB)  for  figure  19,  and  =  1  (SNR  =  0  dB)  for  fig¬ 
ure  20.  For  both  noise  scenarios,  there  is  an  increase  in  AOA  fidelity  (decrease 
in  the  CRLB)  with  increasing  frequency.  This  is  because  the  wavelength  de¬ 
creases  with  increasing  frequency,  thereby  increasing  the  dimensions  of  the 
array  relative  to  a  wavelength.  (Recall  from  sect.  4.1  that  large  arrays,  rela¬ 
tive  to  the  wavelength,  have  greater  phase  differences  between  the  elements, 
so  that  noise  is  mitigated.)  In  the  low  SNR  case,  the  CRLB  is  much  higher, 
particularly  for  short  propagation  distances.  In  fact,  the  CRLB  is  nearly 
independent  of  propagation  distance  for  low  SNR.  This  is  characteristic  of 
a  noise-limited  array,  meaning  that  noise  is  the  dominant  factor  determin- 
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Figure  18:  Sensor  layout  for  circular  RNADS  array  in  a  horizontal  plane. 
Sensors  are  shown  as  dark  circles. 


Figure  19:  CRLB  for  circular  RNADS  array  in  convective  atmospheric 
boundary  layer,  low-noise  case.  SNR  is  20  dB. 

ing  array  performance.  For  high  SNR,  the  array  is  noise  limited  for  short 
propagation  distances,  after  which  the  CRLB  rapidly  increases.  This  rapid 
increase  marks  the  distance  where  the  array  becomes  turbulence  limited.  The 
transition  distance  becomes  shorter  as  the  frequency  is  increased. 

The  near-neutral  atmospheric  case,  with  =  0.01,  is  shown  in  figure  21. 
This  case  looks  nearly  identical  to  the  convective,  high-SNR  case  (fig.  19). 
The  reason  for  this  similarity  is  that  the  values  used  for  in  the  convective 
and  neutral  cases  were  the  same,  and  the  CRLB  responds  primarily  to  small- 
scale  velocity  fluctuations  (the  term  in  eq  (4.5)). 

The  calm  atmospheric  case,  with  =  0.01,  is  shown  in  figure  22.  Al¬ 
though  the  SNR  is  high  in  this  case,  it  looks  qualitatively  similar  to  the 
low-SNR  convective  case  (fig.  20).  This  is  because  both  cases  are  essentially 
noise  limited.  The  turbulence  in  the  calm  atmosphere  is  not  strong  enough 
to  affect  the  CRLB,  at  least  for  the  range  of  firequencies  and  propagation 
distances  shown. 


Figure  20:  CRLB  for  circular  RNADS  array  in  convective  atmospheric 
boundary  layer,  high-noise  case.  SNR  is  0  dB. 


Figure  21;  CRLB  for  circular  RNADS  array  in  neutral  atmospheric  boundary 
layer,  low-noise  case.  SNR  is  20  dB. 


Figure  22:  CRLB  for  circular  RNADS  array  in  calm  atmospheric  boundary 
layer,  low-noise  case.  SNR  is  20  dB. 
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5.  Conclusions 


Turbulence  limits  our  ability  to  obtain  accurate  AOA  estimates  using  acous¬ 
tic  sensor  arrays,  particularly  in  low-noise  environments.  Fortunately,  the 
calculations  in  this  paper  suggest  that  the  AOA  can  be  measured  to  an  ac¬ 
curacy  of  several  degrees  or  better,  even  in  strong  turbulence,  by  the  use  of 
an  array  with  just  a  half  dozen  or  so  sensors. 

Turbulence  effects  on  acoustic  sensor  arrays  are  most  significant  under 
the  following  conditions: 

•  Sufficiently  energetic  turbulence  is  present.  The  type  of  turbulence  that 
predominantly  affects  acoustic  arrays  is  small-scale  turbulence  gener¬ 
ated  by  near-ground  wind  shear  instabilities,  the  strength  of  which  is 
parameterized  by  the  friction  velocity  tt*. 

•  There  is  low  background  noise,  so  that  array  performance  is  limited  by 
turbulence,  rather  than  noise. 

•  Array  dimensions  are  comparable  in  size  to  the  most  energetic  turbu¬ 
lent  eddies.  This  is  usually  true  for  the  arrays  of  interest  in  current 
Army  applications. 

Examples  given  in  this  report  demonstrate  that  the  turbulence  model 
used  to  calculate  array  performance  assessments  can  dramatically  affect  the 
results.  In  particular,  a  von  Karman  turbulence  spectrum  is  generally  prefer¬ 
able  to  a  Gaussian  spectrum,  since  it  models  inertial  subrange  turbulence 
more  realistically. 


26 


6.  Recommendations 


In  this  report,  I  have  only  scratched  the  surface  of  the  many  signal  process¬ 
ing,  turbulence  modeling,  and  acoustic  propagation  issues  relevant  to  design 
and  performance  assessment  of  acoustical  tracking  arrays.  There  are  numer¬ 
ous  areas  for  future  research  that  will  help  improve  our  understanding  and 
enhance  the  performance  of  acoustical  systems: 

1.  Most  importantly,  the  theory  and  predictions  developed  in  this  report 
should  be  tested  by  a  comparison  to  field  data. 

2.  The  propagating  wavefironts  were  assumed  to  be  planar  in  this  report. 
It  would  be  very  desirable  to  extend  the  theory  to  include  spherical 
wavefi-onts,  since  in  many  situations  of  practical  importance  the  sensor 
array  will  be  close  to  the  source. 

3.  The  assumption  of  uncorrelated  noise  at  the  sensor,  as  made  in  sec¬ 
tion  2.2,  is  very  unrealistic.  The  potential  importance  of  using  realistic 
noise  models  needs  to  be  explored. 

4.  Acoustic  tracking  demonstrations  previously  performed  by  ARL  (Srour 
and  Robertson,  1995)  have  often  used  several  separate  independent 
arrays,  each  similar  to  the  one  shown  in  figure  18.  The  main  reason 
for  deploying  multiple  arrays  in  this  fashion  is  to  allow  determination 
of  the  actual  location  (not  just  the  AO  A)  of  the  acoustic  source.  It 
would,  of  course,  be  desirable  to  calculate  performance  estimates  for 
these  multiple-array  deployments.  Performance  assessments  also  are 
needed  when  there  are  multiple  sources. 

5.  Improved  turbulence  models,  including  realistic  inhomogeneity  and 
anisotropy  characteristic  of  atmospheric  turbulence  (Mann,  1994;  Wil¬ 
son,  1996b),  should  be  considered  in  future  research. 

6.  In  field  studies,  acoustic  tracking  has  sometimes  been  found  to  undergo 
periods  of  intermittent  target  loss  (Pham  and  Sadler,  1995).  The  likely 
cause  of  such  dropouts  is  a  failure  of  the  tr2icking  algorithm  brought 
about  by  an  intermittent  change  in  the  meteorological  conditions.  Pre¬ 
viously  (Wilson,  (1996a)),  I  considered  the  effect  of  turbulent  intermit- 
tency  on  acoustic  detection  in  shadows.  Future  research  should  extend 
that  work  to  acoustical  tracking  scenarios. 
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Appendix  A.  Statistical  Turbulence  Models 

In  order  to  determine  the  mutual  coherence  function  (MCF),  and  hence 
model  the  turbulence  effect  on  a  sensor  array,  we  need  a  spectral  model  for 
the  turbulence  field.  The  mathematics  involved  in  properly  defining  and  ma¬ 
nipulating  turbulence  spectra  is  rather  complicated.  The  relevant  equations 
are  summarized  in  this  appendix  and  applied  to  two  different  turbulence 
models:  the  Gaussian  and  von  Karman. 

A-1.  Correlation  and  Spectral  Functions 

Let  us  define  the  correlation  function  of  a  scalar  quantity  c  as 

R(Tiy  T2i  ^3)  =  (^C  (3:^,  X2,  2:3)  C  (2:1,  X2,  3:3)}  1 

where  ri  =  x[—xi,  etc.  For  this  definition,  it  has  been  assumed  that  the  field 
is  homogeneous:  its  statistics  depend  only  on  the  displacement  between  the 
measurement  points.  If  the  field  is  also  isotropic  (i.e.,  its  statistics  are  inde¬ 
pendent  of  coordinate  rotations),  then  the  correlation  function  for  a  scalar 
actually  depends  only  on  the  radial  distance  r  between  the  measurement 
points,  where  =  ri+r^+rl-  In  this  case,  we  can  write  R{r)  =  R  (ri,  r2,  r3) 
for  short. 

The  three-dimensional  spectral  density  function  (spectrum,  for  short)  is 
the  three-dimensional  Fourier  transform  of  the  correlation  function: 

1  TOO  roo  roo 

<^{ki,K2,K3)  =  (A-2) 

exp  [— i  (kiTi  +  K2r2  +  ks^s)]  d^i  dr2  drz, 

/oo  too  too 

/  /  <E>  (ki,  K2,  K3)  (A-3) 

-OO  J  ~oo  J  —  00 

exp  [i  (/ciri  +  K2?’2  +  «3^3)]  dKi  dK2  dwa- 

The  two-dimensional  correlation  function  (TCF,  which  is  required  to  com¬ 
pute  the  MCF  in  sect.  3.2)  is  defined  as  the  two-dimensional  inverse  Fourier 
transform  of  $  (0,  K2,  K3): 

/OO  too 

/  $  (0,  K2,  K3)  exp  [i  {K2r2  +  K3»'3)]  dK2  dK3. 

-00  J -00  (A-4) 

Equivalently,  by  Fourier  transforming  (A-3)  with  respect  to  r\  and  setting 
Ki  =  0,  one  has 

1  t°° 

b(p)  =  :^  /  R{ri,r2,r3)  dri.  (A-5) 

A  note  on  terminology:  Other  authors  have  called  the  TCF  the  transverse 
(rather  than  two-dimensional)  correlation  function.  I  prefer  two-dimensional, 
since  we  derive  the  TCF  by  setting  one  of  the  wavenumbers  in  $  (ki,  K2,  K3) 
to  zero  and  then  taking  the  inverse  Fourier  transform  in  the  remaining  two 
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directions.  Hence  it  is  the  logical  counterpart  of  the  two-dimensional  spec¬ 
trum,  which  we  derive  by  setting  one  of  the  displacements  in  the  correlation 
function  i?  (ri,  r2,  ra)  to  zero,  and  then  taking  the  forward  Fourier  transform 
the  remaining  two  directions. 

Vector  quantities,  such  as  turbulent  velocity  fluctuations,  are  handled 
similarly.  However,  the  situation  becomes  somewhat  more  complicated  for 
vectors,  since  three  directions  are  definable  in  a  given  spectrum  or  corre¬ 
lation:  the  direction  of  the  displacement,  and  also  the  orientations  of  the 
two  velocity  components  under  consideration.  (Only  the  direction  of  the  dis¬ 
placement  is  involved  in  the  scalar  case.)  Let  us  define  the  correlation  as 

Rij  (n,  ^2,  ra)  =  {ui  (x'l,  z'a,  3:3)  uj  {xi,X2,  X3)) ,  (A-6) 

where  the  subscripts  i  and  j  indicate  the  velocity  component.  Note  that 
Rij  is  a  tensor  quantity,  having  nine  components.  By  taking  the  appropriate 
Fourier  transforms  in  a  manner  analogous  to  the  scalar  case,  we  can  define 
tensors  for  the  spectral  density  and  TCF. 


A-2.  Energy  Spectra 

If  only  scalar  quantities  were  of  interest,  it  would  be  a  straightforward  mat¬ 
ter  to  propose  a  model  for  the  correlation  function  or  the  three-dimensional 
spectrum,  and  then  integrate  to  find  the  TCF.  But  in  the  vector  case,  since 
there  are  nine  possible  combinations  for  the  values  of  i  and  j  in  the  correla¬ 
tion  tensor  (eq  (A-6)),  trying  to  model  the  correlation  functions  or  spectral 
density  tensors  directly  can  become  complicated.  It  turns  out  to  be  some¬ 
what  easier  to  start  with  a  model  for  a  scalar  quantity  that  is  a  property  of 
both  scalar  and  vector  fields,  such  as  energy.  For  consistency,  I  adopt  this 
approach  in  both  the  scalar  and  vector  cases.  The  energy  spectrum  E  (k)  is 
defined  such  that  its  integral  over  the  wavenumber  domain  equals  half  the 
total  variance  of  the  field.  For  a  scalar  quantity, 

j^E{K)dK^'^.  (A-7) 

Note  from  the  Fourier  transform  relation  (eq  (A-3)),  with  n  =  r2  =  rs  =  0, 

/oo  roo  poo 

/  /  $  («i,  K2^  «3)  d«2  dK3. 

-OO  J  —  oo  J  —  oo 

Hence,  we  find  the  three-dimensional  spectrum  $  (k)  by  multiplying  E  («) 
by  two,  and  then  dividing  by  the  “area”  of  a  spherical  shell  in  wavenumber 
space, 

=  (A.8) 

Substituting  into  the  definition  of  the  TCF  (eq  (A-4)),  we  obtain 


exp  [i  (^2^2  +  K'Zf’z)]  dK2  dKz, 

(A-9) 


where  +  rl- 

For  vectors,  such  as  the  turbulent  velocities,  E  (k)  is  defined  such  that 


(A-10) 


If  this  definition  is  used,  represents  the  variance  in  just  one  of  the  three 
velocity  components.  The  three-dimensional  spectra  are  related  to  the  energy 
spectrum  according  to  (Batchelor,  1953) 


From  the  above,  the  reader  can  verify,  analogously  to  the  scalar  case,  that 
$11  (k)  -t-  $22  (k)  +  ^33  («)  = 


Substituting  $ii  (k),  as  given  by  equation  (A-11),  into  the  definition  of  the 
TCF  (eq  (A-4)),  we  have 


1  TOO  fOi 


(y/ 

K2  +  «3 


exp  [i  (k2^'2  +  «3J’3)]  <^«2  dK3. 

(A-12) 


Quite  interestingly,  this  result  is  exactly  half  the  result  for  a  scalar.  Hence, 
if  we  choose  the  same  functional  form  for  the  energy  spectrum  in  both  the 
scalar  and  vector  cases,  the  resulting  TCFs  will  likewise  have  the  same  func¬ 
tional  form. 


A-3.  Length  Scales 

It  is  important  that  we  be  able  to  accurately  quantify  length  scales  associated 
with  given  spectra.  By  length  scale,  I  mean  a  value  representative  of  the  size 
of  the  most  energetic  eddies  in  the  spectrum.  For  a  scalar  quantity,  the 
integral  length  scale  is  defined  as 

£  =  “2  /  R{r)  dr.  (A-13) 

(T‘‘  Jo 

It  turns  out  that  the  integral  length  scale  and  TCF  are  simply  related.  From 
the  Fourier  transform  relation  between  the  correlation  function  and  the  spec¬ 
trum,  and  from  equation  (A-8),  we  have 


Setting  r2  =  rs  =  0,  k\  =  0,  and  using  the  fact  that  J?(ri,0,0)  is  an  even 
function,  we  obtain 


1. 

TT 


foo  1  roo  roc 


dK2  dK3. 
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(A-14) 


Comparing  now  with  equation  (A-9),  we  find  that 

/:  =  ^6(o). 

In  the  vector  case,  the  direction  of  the  velocity  components  as  well  as  the 
direction  of  the  displacement  must  be  considered.  Two  significant  scales  can 
be  defined:  one  where  the  displacement  is  parallel  to  the  velocity,  e.g., 

1  r°° 

^\\  =  ~2  Rn{ri,0,0)dri, 
cr  Jo 

and  the  other  where  it  is  perpendicular,  e.g., 

1  r°° 

£x  =  “2  /  ■R22(’"i,0,0)  dri. 

<7"“  Jo 

For  homogeneous,  isotropic  turbulence,  it  can  be  proven  that  (Batchelor, 
1953) 


£||  —  2/lx- 


(A-15) 


Furthermore,  using  a  derivation  similar  to  the  scalar  case,  one  finds 

£||  =  ^bn  (0) .  (A-16) 

IF 


4.  Gaussian  Model 

Let  US  now  consider  specific  models  for  atmospheric  turbulence.  Unfortu¬ 
nately,  we  do  not  now  have  (and  probably  never  will  have)  three-dimensional 
turbulence  models  that  are  known  to  work  satisfactorily  for  a  variety  of  at¬ 
mospheric  conditions.  Hence,  we  are  forced  to  consider  idealized  models. 
In  this  section,  I  describe  one  such  idealization,  the  Gaussian  model.  The 
main  advantage  of  the  Gaussian  model  is  its  analytical  convenience;  unfor¬ 
tunately,  it  can  work  well  for  qualitatively  modeling  only  the  largest  scale 
features  of  the  turbulence  (called  the  energy-containing  subrange).  The  von 
Karman  model,  described  in  section  A-5,  is  somewhat  more  general,  since  it 
also  describes  well  the  smaller  scale  eddies  (the  inertial  subrange). 

It  is  customary  to  define  the  Gaussian  model  in  terms  of  its  three- 
dimensional  correlation  function  R{r).  For  the  scalar  case, 

R{r)  =  cr^exp  >  (A-17) 

where  L  is  a  length  scale.  The  corresponding  three-dimensional  spectrum, 
found  by  Fourier  transformation  of  the  correlation  function  in  all  three  di¬ 
mensions,  is 

,  .  (  k^L^\  /a 

=  j-  (A-18) 

The  energy  spectrum,  which  follows  from  multiplication  by  27r«;^,  is  shown 
in  figure  A-1. 


Figure  A-1:  Comparison  of  Gaussian  and  von  Karman  energy  spectral  mod¬ 
els. 


The  TCF  is  found  from  the  energy  spectrum  and  equation  (A-9).  The 
integration  is  not  difficult  for  the  Gaussian  function,  with  the  result  being 

Hence  from  equation  (A-14),  we  have 

£  =  ^L.  (A-20) 

We  find  the  TCF  for  the  velocity  field  by  multiplying  by  3  (to  account  for 
the  vector  nature  of  this  quantity,  see  eq  (A-10)),  and  then  dividing  by  2, 
as  discussed  earlier  in  connection  with  equation  (A-12).  By  equation  (A-16), 
then. 


(A-21) 


A-5.  Von  Karman  Model 


The  von  Karman  type  of  model  is  developed  directly  from  the  following 
proposed  form  for  the  energy  spectrum  of  a  scalar: 


4r(i/ -1-5/2) 

Z^r{u)  (1  +  «2^2)«'+5/2- 


(A-22) 


The  parameter  £  is  a  characteristic  length  scale,  and  i/  controls  the  power- 
law  dependence  in  the  inertial  subrange  (k£  :§>  1).  Generally,  we  set  u  = 
1  /3  to  obtain  Kolmogorov’s  k"®/®  power  law  for  the  inertial  subrange.  The 
von  Karman  model  is  compared  to  the  Gaussian  model  in  figure  A-1.  The 
main  difference  between  the  von  Karman  and  Gaussian  models  is  that  the 
Gaussian  model  decays  much  more  rapidly  at  large  wavenumbers. 
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Calculating  the  TCF  from  relationship  (A-9),  we  find 


b{p)  = 


23/2-t/^2^ 

(i/) 


1/+1/2  r 


^v+1/2 


(A-23) 


Integrals  (3.773.6)  and  (6.726.4)  in  Gradshteyn  and  Rhyzhik  (1994)  were 
used  to  derive  this  result.  Ki^  is  the  modified  Bessel  function  of  the  second 
kind.  As  discussed  above,  the  result  bn  (p)  for  vectors  is  simply  3/2  times 
the  scalar  result. 

Integral  length  scales  can  be  derived  by  replacement  of  the  Bessel  func¬ 
tions  by  the  first  term  in  their  Taylor  series  expansions.  The  result  for  scalars 
is 


C  = 


2^T{u  +  l/2) 

3r(i.) 


For  vectors, 


v^r  (i/  + 1/2) 

Hi  =  — rr:A — 


(A-24) 


(A-25) 


Prom  equation  (A-20)  and  (A-24),  it  can  be  shown  that  L  for  the  Gaus¬ 
sian  model  is  1.78  times  i  for  the  von  Karman  model,  assuming  the  same 
value  for  L  is  used  in  each  case. 
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